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First, let's go over the important steps in the process of writing a SPICE-based 
Fortran program and putting it to work: 


¢ Understand the geometry problem. 


¢ Identify the set of SPICE kernels that contain the data needed to perform the 
computation. 


¢ Select the SPICE APIs needed to compute the quantities of interest. 
¢ Write and compile the program. 


¢ Get actual kernel files and verify that they contain the data needed to support 
the computation for the time(s) of interest. 


¢ Run the program. 


To illustrate these steps, let's write a program that computes the apparent 
intersection of the boresight ray of a given CASSINI science instrument with the 
surface of a given Saturnian satellite. The program will compute 


¢ Planetocentric and planetodetic (geodetic) latitudes and longitudes of the 
intercept point. 


¢ Range from spacecraft to intercept point. 
¢ Illumination angles (phase, solar incidence, and emission) at the intercept point. 
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Which Kernels are Needed? 
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Data required to compute vectors, rotations and other parameters shown in 
the picture are stored in the SPICE kernels listed below. 


Note: these kernels have been selected to support this presentation; they should not be assumed to be 


appropriate for user applications. 


Parameter 


time conversions 


satellite orientation 
satellite shape 
satellite position 


planet barycenter position 
spacecraft position 
spacecraft orientation 
instrument alignment 
instrument boresight 
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Kernel Type 
generic LSK 
CASSINI SCLK 
CASSINI PCK 
CASSINI PCK 
planet/sat 
ephemeris SPK 
planet SPK 
spacecraft SPK 
spacecraft CK 
CASSINI FK 
Instrument IK 
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The easiest and most flexible way to make required kernels available to the 
program is via FURNSH. For this example we make a setup file (also called a 
“metakernel” or “furnsh kernel’) containing a list of kernels to be loaded: 


Note: these kernels have been selected to support this presentation they 
should not be assumed to be appropriate for user applications. 
\begindata 


KERNELS TO LOAD = ( 


\begintext 


and we make the program prompt for the name of this setup file: 


CALL PROMPT ( 'Enter setup file name > ', SETUPF ) 
CALL FURNSH ( SETUPF ) 
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¢ Prompt for setup file (“metakernel”) name load kernels specified via setup 
file. (Done on previous chart.) 


¢ Prompt for user inputs required to completely specify problem. Obtain 
further inputs required by geometry routines via SPICELIB calls. 


¢ Compute the intersection of the boresight direction ray with the surface of 
the satellite, presented as a triaxial ellipsoid. 


If there is an intersection, 
-Convert Cartesian coordinates of the intercept point to planetocentric 
latitudinal and planetodetic coordinates 
«Compute spacecraft-to-intercept point range 
¢Find the illumination angles (phase, solar incidence, and emission) at 
the intercept point 


¢ Display the results. 


We discuss the geometric portion of the problem next. 
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Compute the intercept point (Pornt) of the boresight vector (INsITE) specified in 
the instrument frame (IFRAME) of the instrument mounted on the spacecraft (scnm) 
with the surface of the satellite (satnm) at the TDB time of interest (ET) in the 
satellite’s body-fixed frame (FIxREF). This call also returns the light-time 
corrected epoch at the intercept point (TRcEPc), the spacecraft-to-intercept point 
vector (SRFVEC), and a flag indicating whether the intercept was found (Founp). 

We use "converged Newtonian" light time plus stellar aberration corrections to 
produce the most accurate surface intercept solution possible. We model the 
surface of the satellite as an ellipsoid. 


CALL SINCPT ( 'Ellipsoid', SATNM, ET, FIXREF, 'CN+S', SCNM, IFRAME, 
INSITE, POINT, TRGEPC, SRFVEC, FOUND ) 


The range we want is obtained from the outputs of stncpt. These 

outputs are defined only if a surface intercept is found. If Founp is true, the 
spacecraft-to-surface intercept range is the norm of the output argument SRFVEc. 
Units are km. We use the SPICELIB function vnors to obtain the norm: 


VNORM ( SRFVEC ) 


We'll write out the range data along with the other program results. 
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Compute Lat/Lon and Illumination Angles 
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Compute the planetocentric latitude (PcLAT) and longitude (PcLon), as well as 
the planetodetic latitude (PpLAT) and longitude (PpLoN) of the intersection 


point. 
IF ( FOUND ) THEN 
CALL RECLAT ( POINT, R, PCLON, PCLAT ) 
Cc Let RE, RP, and F be the satellite's longer equatorial 
C radius, polar radius, and flattening factor. 
RE = RADII(1) 
RP = RADII (3) 


CALL RECGEO ( POINT, RE, F, PDLON, PDLAT, ALT ) 
The illumination angles we want are the outputs of ILumin. Units are radians. 


CALL ILUMIN ( 'Ellipsoid', SATNM, ET, FIXREF, 
"CN+S', SCNM, POINT, TRGEPC, SRFVEC, 


PHASE, SOLAR, EMISSN ) 
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Geometry Calculations: Summary 
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Cc Compute the boresight ray intersection with the surface of the 
Cc target body. 


CALL SINCPT ( 'Ellipsoid', SATNM, ET, FIXREF, 'CN+S', SCNM, IFRAME, 
INSITE, POINT, TRGEPC, SRFVEC, FOUND ) 


Cc If an intercept is found, compute planetocentric and planetodetic 
Cc latitude and longitude of the point. 


IF( FOUND ) THEN 


CALL RECLAT ( POINT, R, PCLON, PCLAT ) 


Cc Let RE, RP, and F be the satellite's longer equatorial 
Cc radius, polar radius, and flattening factor. 

RE = RADII(1) 

RP = RADII(3) 

F = (RE - RP) / RE 


CALL RECGEO ( POINT, RE, F, PDLON, PDLAT, ALT ) 
Cc Compute illumination angles at the surface point. 


CALL ILUMIN ( 'Ellipsoid', SATNM, ET, FIXREF, 'CN+S', SCNM, 
POINT, TRGEPC, SRFVEC, PHASE, SOLAR, EMISSN ) 


ELSE 
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Get Inputs - 1 
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The code above used quite a few inputs that we don't have yet: 


¢ TDB epoch of interest (ET ) 
¢ satellite and s/c names (SATNM, SCNM) 


¢ satellite body-fixed frame name (FrxR=EF) 


¢ satellite ellipsoid radii (RADII) 
¢ instrument fixed frame name (IFRAME) 
¢ instrument boresight vector in the instrument frame (INSITE) 


Some of these values are user inputs others can be obtained via SPICELIB calls 


once the required kernels have been loaded. 


Let's prompt for the satellite name (satny), satellite frame name (FIXREF), 
spacecraft name (scnwy), instrument name (INstTNM) and time of interest (TIME): 


CALL 
CALL 
CALL 
CALL 
CALL 


PROMPT 
PROMPT 
PROMPT 
PROMPT 
PROMPT 
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Then we can get the rest of the inputs from SPICELIB calls: 


To get the TDB epoch (ET) from the user-supplied time string (which may 
refer to the UTC, TDB or TT time systems): 


CALL STR2ET ( TIME, ET ) 
To get the satellite’s ellipsoid radii (RADII): 


CALL BODVRD ( SATNM, 'RADII', 3, I, RADII ) 


To get the instrument boresight direction (INSITE) and the name of the 
instrument frame (IFRAME) in which it is defined: 


CALL GETFVN ( INSTNM, ROOM, SHAPE, IFRAME, 
INSITE, N, BUNDRY ) 
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Getting Inputs: Summary 
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Cc Prompt for the user-supplied inputs for our program. 
CALL PROMPT ( 'Enter setup file name > ', SETUPF ) 
CALL FURNSH ( SETUPF ) 


CALL PROMPT( 'Enter satellite name > ', SATNM ) 

CALL PROMPT( 'Enter satellite frame > ', FIXREF ) 

CALL PROMPT( 'Enter spacecraft name > ', SCNM ) 

CALL PROMPT( 'Enter instrument name > ', INSTNM ) 

CALL PROMPT( 'Enter time > ', TIME ) 
Cc Get the epoch corresponding to the input time: 


CALL STR2ET ( TIME, ET ) 


Cc Get the radii of the satellite. 
CALL BODVRD ( SATNM, 'RADII', 3, I, RADII ) 
Cc Get the instrument boresight and frame name. 
CALL GETFVN ( INSTNM, ROOM, SHAPE, IFRAME, 
INSITE, N, BUNDRY ) 
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Display Results 
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Cc Display results. Convert angles from radians to degrees 

Cc for output. 
WRITE ( *, '(1X,A,F12.6)' ) 
"Intercept planetocentric longitude (deg) : ',  DPR() *PCLON 
WRITE ( *, '(1X,A,F12.6)' ) 
"Intercept planetocentric latitude (deg) : ', DPR()*PCLAT 
WRITE ( *, '(1X,A,F12.6)' ) 
"Intercept planetodetic longitude (deg) : ',  DPR() *PDLON 
WRITE ( *, '(1X,A,F12.6)' ) 
"Intercept planetodetic latitude (deg) : ', DPR()*PDLAT 


WRITE ( *, '(1X,A,F12.6)' ) 

"Range from spacecraft to intercept point (km): un 

VNORM (SREVEC) 

WRITE ( *, '(1X,A,F12.6)' ) 

"Intercept phase angle (deg) : ',  DPR() *PHASE 
WRITE ( *, '(1X,A,F12.6)' ) 

"Intercept solar incidence angle (deg) : ', DPR() *SOLAR 
WRITE ( *, '(1X,A,F12.6)' ) 

"Intercept emission angle (deg): ', 

DPR () *EMISSN 


ELSE 
WRITE (*,*) 'No intercept point found at '// TIME 
END IF 
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To finish up the program we need to declare the variables we've used. 


¢ We'll highlight techniques used by NAIF programmers 
¢ Add remaining Fortran code required to make a syntactically valid program 
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PROGRAM PROG26 
IMPLICIT NONE 


DOUBLE PRECISION 
DOUBLE PRECISION 


INTEGER 
PARAMETER 
INTEGER 
PARAMETER 
INTEGER 
PARAMETER 


CHARACTER* (WORDSZ) 
CHARACTER* (WORDSZ) 
CHARACTER* (WORDSZ) 
CHARACTER* (WORDSZ) 
CHARACTER* (WORDSZ) 
CHARACTER* (FILESZ) 
CHARACTER* (WORDSZ) 
CHARACTER* (WORDSZ) 


DOUBLE PRECISION 
DOUBLE PRECISION 
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Complete Source Code - 1 
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DPR 
VNORM 


FILESZ 

( FILESZ 
WORDSZ 

( WORDSZ 
ROOM 

( ROOM 


IFRAME 
INSTNM 
SATNM 
FIXREF 
SCNM 
SETUPF 
SHAPE 
TIME 


ALT 


255 


40 
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BUNDRY (3, ROOM) 


) 


) 
) 


DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 


PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 


INTEGER 
INTEGER 


LOGICAL 


EMISSN 
ET 

F 

INSITE (3) 
SRFVEC (3) 
PCLAT 
PCLON 
PDLAT 
PDLON 
PHASE 
POINT (3) 
R 

RADII (3) 
RE 

RP 

SOLAR 
TRGEPC 


FOUND 
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Cc Prompt for the user-supplied inputs for our program. 
CALL PROMPT ( 'Enter setup file name > ', SETUPF ) 
CALL FURNSH ( SETUPF ) 
CALL PROMPT ( 'Enter satellite name > ', SATNM ) 
CALL PROMPT ( 'Enter satellite frame > ', FIXREF ) 
CALL PROMPT ( 'Enter spacecraft name > ', SCNM ) 
CALL PROMPT ( 'Enter instrument name > ', INSTNM ) 
CALL PROMPT ( ‘Enter time > ', TIME ) 

Cc Get the epoch corresponding to the input time: 
CALL STR2ET ( TIME, ET ) 

C Get the radii of the satellite. 
CALL BODVRD ( SATNM, 'RADII', 3, I, RADII ) 

Cc Get the instrument boresight and frame name. 


CALL GETFVN ( INSTNM, ROOM, SHAPE, 
INSITE, 
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Complete Source Code - 2 
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N, BUNDRY 


IFRAME, 
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Complete Source Code - 3 
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Compute the boresight ray intersection with the surface of the 
target body. 
CALL SINCPT ( 'Ellipsoid', SATNM, ET, FIXREF, 'CN+S', SCNM, IFRAME, 
INSITE, POINT, TRGEPC, SRFVEC, FOUND ) 

If an intercept is found, compute planetocentric and planetodetic 
latitude and longitude of the point. 
IF( FOUND ) THEN 

CALL RECLAT ( POINT, R, PCLON, PCLAT ) 

Let RE, RP, and F be the satellite's longer equatorial 

radius, polar radius, and flattening factor. 


RE = RADII(1) 
RP = RADII(3) 


CALL RECGEO ( POINT, RE, F, PDLON, PDLAT, ALT ) 
Compute illumination angles at the surface point. 


CALL ILUMIN ( 'Ellipsoid', SATNM, ET, FIXREF, 'CN+S', SCNM, 
POINT, TRGEPC, SRFVEC, PHASE, SOLAR, EMISSN ) 

Display results. Convert angles from radians to degrees 

for output. 

WRITE ( *, * ) 

WRITE ( *, '(1X,A,F12.6)' ) 

"Intercept planetocentric longitude (deg) : ',  DPR() *PCLON 
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Complete Source Code - 4 
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WRITE ( *, '(1X,A,F12.6)' ) 


‘Intercept planetocentric latitude (deg) : ', DPR()*PCLAT 
WRITE ( *, '(1X,A,F12.6)' ) 
"Intercept planetodetic longitude (deg) : ',  DPR() *PDLON 
WRITE ( *, '(1X,A,F12.6)' ) 
"Intercept planetodetic latitude (deg) : ', DPR()*PDLAT 


WRITE ( *, '(1X,A,F12.6)' ) 

"Range from spacecraft to intercept point (km): ', 

VNORM (SRFVEC) 

WRITE ( *, '(1X,A,F12.6)' ) 

"Intercept phase angle (deg) : ', DPR() * PHASE 
WRITE ( *, '(1X,A,F12.6)' ) 

"Intercept solar incidence angle (deg) : ',  DPR() *SOLAR 
WRITE ( *, '(1X,A,F12.6)' ) 

"Intercept emission angle (deg): ', 

DPR () *EMISSN 


ELSE 

WRITE (*,*) 'No intercept point found at '// TIME 
END IF 
END 
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¢ First be sure that both the SPICE Toolkit anda 
Fortran compiler are properly installed. 


— A “hello world" Fortran program must be able to compile, link, 
and run successfully in your environment. 


— Any of the mkprodct.* scripts in the toolkit/src/* paths of the 
SPICE Toolkit installation should execute properly. 


¢ Ways to compile and link the program: 


— If you're familiar with the "make" utility, create a makefile. Use 
compiler and linker options from the mkprodct.* script found in 


the toolkit/src/cookbook path of your SPICE Toolkit installation. 


— Or, modify the cookbook mkprodct.* build script. 


» Your program name must be *.pgm, for example demo.pgm, 
to be recognized by the script. 


» Change the library references in the script to use absolute 
pathnames. 


» Change the path for the executable to the current working 
directory. 


» On some platforms, you must modify the script to refer to 
your program by name. 
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Compile and Link the Program - 2 
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— Or, compile the program on the command line. The program 
must be linked against the SPICELIB object library spicelib.a 
(spicelib.lib under MS Windows systems). Ona PC running 
Linux and g77, if 


» The g77 compiler is in your path 
¢ As indicated by the response to the command "which g77" 


» the Toolkit is installed in the path (for the purpose of this 
example) /myhome/toolkit 


» You've named the program demo.f 


then you can compile and link your program using the 
command 


» g77 -o demo demo.f \ 
/myhome/toolkit/lib/spicelib.a 
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Compile and Link the Program - 3 
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Terminal Window 


Prompt> mkprodct..csh 


Using the g7/77 compiler. 


Setting default Fortran compile options: 
=. = 


Setting default C compile options: 
=0 


Setting default link options: 


Compiling and linking: demo.pgm 
Compiling and linking: demo.pgm 


PrOmp TC 
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Running the Program - 1 
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It looks like we have everything taken care of: 


¢ We have all necessary kernels 


¢ We made a setup file (metakernel) pointing to them 


¢ We wrote the program 


¢ We compiled and linked it 


Let's run it. 
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Running the Program - 2 
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Terminal Window 


Prompt> demo 
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Backup 
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¢ Latitude definitions: 


— Planetocentric latitude of a point P: angle between segment from 
origin to point and x-y plane (red arc in diagram). 


— Planetodetic latitude of a point P: angle between x-y plane and 
extension of ellipsoid normal vector N that connects x-y plane and 
P (blue arc in diagram). 


z-axis 
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